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ABSTRACT 


In this thesis, we will study bio-mathematics. We will introduce differential equations, 
biological applications, and simulations with emphasis in molecular events. One of 
the first courses of action is to introduce and construct a mathematical model of our 
biological element. The biological element of study is the Hepatitis C virus. The idea 
in creating a mathematical model is to approach the biological element in small steps. 
We will first introduce a block (schematic) diagram of the element, create differential 
equations that define the diagram, convert the dimensional equations to non-dimensional 
equations, reduce the number of parameters, identify the important parameters, and 
analyze the results. These results will tell us which variables must be adjusted to 


prevent the Hepatitis C virus from becoming chronic. 
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Introduction 


Mathematical models are useful when there is reason to believe that emergent 
properties constitute the system. These qualitative features are earmarked when policy 
changes must be reached. The decision to change policy may involve both qualitative 
and quantitative analysis that are used to predict behavior under certain conditions or 
decide which parameters enhance the spread of disease, and hence, which public actions 
should be taken to counter the effect of that disease. These qualitative features may also 
be used to calculate the number of vaccines required to eradicate the disease or, at least, 
get it under control. An example of this type of mathematical model is presented in 
the article by Marco Arieli Herrera-Valdez, Maytee Cruz-Aponte and Carlos Castillo- 
Chavez entitled, “Multiple Outbreaks for the Same Pandemic: Local Transportation 
and Social Distancing Explain the Different “Waves” of A-HIN1IPDM Cases Observed 
in Mexico during 2009”. As we look closely at the steps to formulate a mathematical 
model, based on the material from the book, ”A Primer on Mathematical Models in 
Biology” written by Lee A. Segel and Leah Edelstein-Keshet, we may reference this 


article. 


Chapter 1 


Background of a Mathematical 
Model in Biology 


1.1 Introduction of a Model 


As Segal and Edelstein-Keshet (2013) states, ”A model can be described as 
a caricature of a real system”. A good caricature captures the essence of the system 
and neglects the detail of the system. Establishing this essence from the complexity 
of the system is the key to forming an informative mathematical model. A descriptive 
mathematical model is the simplest type, such as experimental observations that can 
be approximated to a straight line or a sum of exponential functions. A computer sim- 
ulation is another type of model that increases the level of detail and brings faithfulness 
to the underlying model. There are three distinct reasons for using simulations. Firstly, 
simulations help the investigator get an initial “feel” for the behavior of a model be- 
fore spending valuable time in the details of the analysis. Secondly, simulations help 
the investigator make accurate quantitative predictions. Thirdly, simulations help the 
investigator explore advanced versions of models that are not very easy to analytically 
track. The ability to pick a model type is an important step. A model can be viewed 
as a lie since details can be neglected and possible important features distorted to give 
rise to the essential aspects of the model. We must keep in mind that an initially wrong 
model should not be rejected, and an initially right model should not be accepted. As 


Picasso says of art, we must think of in a model, ”a lie that helps us see the truth”. 


[SEK 13] 


1.1.1 Schematic Diagram 


In model building, a crucial first step is conception of a verbal model or a 
schematic diagram. The thoughts describing the underlying mechanism and/or the 
words used to describe the relationship between parts of a system maintain a focus 
on the choice of experiments to be used. From this step, formulating a mathematical 
model is easier. The best process for this step is to amass an accurate depiction of the 
relevant background. We need to ask ourselves what the essential ingredients in the em- 
bodiment of the phenomenon under investigation are. We need to list the “unknowns” 
or dependent variables and their relationship to one another. The main independent 
variable in almost all biological investigations is time. The next step is to determine if 
the phenomenon is regarded as probabilisitic or deterministic. If deterministic, we use 
differential equations and regard continuous changes as discrete. The hardest step in 
formulating a mathematical model is to produce a proper set of equations. The most 
important aid for this step is bookkeeping. The process of bookkeeping keeps track 
of some quantity while others remain invariant. In biological problems, this process is 


difficult and often requires more difficult steps to be taken. [SEK13] 


1.1.2 Solving Formulated Equations 


After equations are formulated, we must solve them. Most often a combina- 
tion of numeric, geometric, and analytic approaches are utilized for this task. Software 
programs such as MAPLE, MATLAB, Mathematica, and XPP are used to approximate 
numerical solutions for a number of representative parameter values. These numerical 
solutions can be used in preliminary simulations to show that a developed theory can 
explain appropriate sets of data. Analytical results are often required when trying to 
determine if the represented parameter sets agree with the desired result or if they are 


made to agree with the desired result. [SEK13] 


1.1.3. Qualitative Results 


Quantitative and qualitative results are two types of conclusions found when 
the model equations are manipulated. The qualitative conclusions are useful in the 
identification of solutions after a long period of time, called attractors. The attractor 
solutions are typically a steady state or periodic oscillation over time. These types of 
solutions often lead to a “rule of thumb” such as the time affiliated with a chemical 
reaction, t + 1/k, where k is the rate constant at which a substance A breaks down to 


substance B. [SEK13] 


It is possible to find qualitative conclusions without prior knowledge of the 
magnitude of present parameters. However, when models start to become complex it is 
useful to obtain rough magnitude estimates for those various parameters. Such estimates 
can be researched from prior literature or developed from general intuition. The most 
desirable way to obtain these estimates is when the values of all parameters can be 
determined from experiments other than the current experiment in question. Then the 
current experiment can be used to verify or disprove the validity of our model. Altering a 
large number of parameters to fit a variety of different experimental results is difficult. 
This is true because models are, typically, nonlinear and changing one parameter to 


“fix-up” a deficiency can make the previously obtained results inconclusive. [SEK13] 


1.1.4 Robustness and Analysis 


Although the development of specific experimental predictions is important, 
the concepts that the conclusions produce are advantageous in designing new experi- 
ments for further analysis. Therefore, the robustness of the model can be even more 
important than the agreement with the specific experiment. The principal conclusions 


of the model must be maintained whether details of the model are amended or not. 


Once a model is constructed, one must analyze the results and find the key 
features hidden that yield the major conclusions. These conclusions may not be recep- 
tive to experimental tests, but are still valuable in developing concepts. Models should 


be considered guides to acquiring experimental information. For a complex system, no 


single model is appropriate and different questions require different models. [SEK13] 


1.2. Biochemical Kinetics 


In this next section, we will write the governing differential equations of bio- 
chemical kinetics by considering several examples of chemical dynamics and presenting 
several principles that can be used to simplify these equations. After these principles 
are applied, we will see the utility of these simplifications and show how simulations can 


complement the analysis of a mathematical model. 


1.2.1 Kinetic Scheme 


Before we begin, let us take a look at the meaning of a kinetic scheme. A 
kinetic scheme is the conventional description of a molecule shifting between two states, 
A and B, where A and B represent the concentrations (number per unit volume) of two 


molecular configurations. Shown as 


An ete. - (oT) 


where k, and k_, are the rate coefficients. The rate constant, k;, times a small time 
interval, At, is defined as the probability that a molecule in state A shifts to state 
B and stays in state B. Similarly, if At is sufficiently small, then k_,At is a good 
approximation to the probability that a molecule initially in state B changes to state 
A and remains in state A. We observe that our definitions imply that k,,k_, have 
units of 1/time and thus, their reciprocals are characteristic time scales. The time step 
should be smaller than the characteristic time for the transition between states. Then 
At < 1/k; and At < 1/k_1. We will follow the Markov properties for the shift of a 
molecule between two configurations. They are: 

(M1) Transitions between states are random. 

(M2) The probability that a transition occurs during some time interval does not de- 


pend on the history of events preceding the time in question. 


(M3) If environmental conditions are fixed then the overall characteristics of the tran- 
sitions that occur in some time interval, do not depend on the time at which the obser- 


vations are made. [SEK13] 


1.2.2 Differential Equations 


Now, let us derive differential equations for the change in time of concentrations 
A and B. If there are A molecules per unit volume, the expected decrease in the number 
of these molecules during At is: 
decrease in A molecules = total number of A molecules x fraction that becomes B 
= A x (k, At) 


Thus the following equation describes the expected change in the number of 


A molecules during the time interval (t,t + At): 


A(t + At) — A(t) = —A(t) - ky At + B(t) - (kK) At. 


Dividing by At, taking the limit as At > 0, and using the definition of derivatives leads 


to: 


A B 
(2.2a) a =-k,jA+k_1B and “ =k,A—k_1B. (2.26) 


This gives us two differential equations describing the kinetic scheme of A 
and B. However, before we move forward, we need to complete the mathematical 
translation by prescribing the initial state of the system at time t = 0. Let A(0) = Ao 
and B(0) = Bo, then the differential equations and the initial conditions make up the 
model. We can simplify this model by using the conservation law, A(t) + B(t) = M, 
noting that M is a constant since the molecules only shift between conformations and 
do not degrade over such a short period of time. Thus, M = Ag+ Bo, the total number 
of molecules at t = 0. Using this conservation law at time t we can solve for B(t) and 


substitute into (2.2a). In doing so, we get 


= —kA+h(M — A), 


which leads to 


dA 

i —(ky + k_1)A +kiM. 

dt 

Now, solving for A requires a review of first order differential equations. Con- 
sider a more general version of an = —k,A, such as 

dx 
—=k(t 
T= MO 


with initial condition x(0) = xo. This equation is called a First-Order Differential 
Equation, ODE, since the highest derivative in the equation is the first derivative. The 
general solution is 


t 
x(t) = Cexp(K(t))t, | where C is a constant and K(t) = | k(s)ds. 
0 


This solution is derived as follows: 


Upon integrating over the interval 0 < x < t, xo < x(s) < x(t), we obtain: 


t 
In(u)|2{¢) = In(a(t)) — In(w()) = In() — In(a) = In(=) = | k(s)ds. 
0 0 

Exponentiation of both sides yields the solution 

t 

x(t) = xpexp / k(s)as]. 
0 

In the case that k(t) is constant, say k(t) = r, we evaluate the integral and obtain 


ere | / | Sone Raa 


Therefore, it can be shown that all solutions of ae = k(t)x are special cases of 
a(t) = xoexpK (t)t for certain values of ao. Thus, the solution for 4 = —(ky +k-1)A+ 


k_1M with initial condition Ao is 


k_1M 
—_ kiM 
In the above, C= Ao = Aso and Ass = ER 


Substituting, we have 


A(t) = Ago — (Aco — Ape thE, 


Similarly, 


B(t) = M — A(t) = (M — Ax) + (Aco — Ao)e7 tt, 


Notice as t + oo, A(t) approaches wots and B(t) approaches CEEEE 


At steady state, 


dB 
—=—k,A {B= —=kA-k,jB= 
dt af + ky 0 and dt 1 k_y 0, 


which occurs when the concentrations of A and B are constant. Here the rate of con- 
version of A to B should exactly balance the rate of conversion of B to A. Thus, 
kj A = k_,1B. The time it takes for A and B to reach their steady states is the time 
scale (k, + k_1)~'. [SEK13] 


1.2.3. Deterministic versus Stochastic Approaches to Solving Problems 


Now, we will consider deterministic versus stochastic approaches to solving 
problems related to biochemical kinetics. A basic assumption of our model is a proba- 
bilistic transition between the two states. Using Newton’s laws one can derive a deter- 
ministic problem, in which, its solution describes the gross motion of molecules. Models 
of this kind have already been developed and form a basis for numerical simulations 


that give information about large molecule dynamics. [SEK13] 


Our equation (2.1) describing the expected change in the number of A molecules 


during time interval (t + At) is: 


A(t + At) — A(t) = —A(t)- (kr, At) + B(t)- (kA), — (2.3) 


which describes probabilistic assumptions for the change in the expected or average 
number of A molecules. Monte Carlo simulations can be used to depict the stochastic 
shift between A and B states. Suppose that at time, t, a molecule is in state A. We 
know that there is a probability, k,; At, that during the time interval (t,t + At) it will 
shift to the B state. If we select a random number between 0 and 1 and that number 
is between 0 and k, At then the simulation shifts the configuration to B; otherwise the 
configuration remains A. Repeating this calculation, one can develop a simulated his- 
tory of the states, A and B, of a single molecule during some time interval. This history 
can be used to compare simulated data with experiments. There are, also, analytical 
methods that will allow conclusions to be drawn if one retains the stochastic character 


of (2.1). [SEK13] 


We have considered a type of phenomena where an average number of molecules 
shift from one state to another. Let’s consider a phenomena where we consider a shift in 
a single channel molecule from an open state to a closed state. These channel molecules 
are responsible for the electrical conductance through cell membranes. Model (2.3) re- 
mains relevant for such situations if A and B are interpreted as the probability that 
a single channel is, respectively, open and closed. Single channel recordings yield in- 
formation not only on mean values of A(t) and B(t), but also on standard deviations 
and other statistical measures. The following graph, Figure 1.1, shows observation of a 


single channel shifting from its open to closed states. [SEK13} 
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Figure 1.1: Channel Shifting from Open to Closed States. [SEK13] 


The derivative of A(t) in this case is either zero or infinite. This charac- 
terization of the derivative remains true no matter how many channels are observed. 
However, if many channels are monitored the jumps between states become less notice- 
able, as shown in the following graph, Figure 1.2, then a true jumpy curve can be well 
approximated by a smooth curve. This smooth curve will have a well-behaved deriva- 
tive, and it is this curve that we seek when we solve for A(t) in our differential equation 


formulation of kinetics. [SEK13] 


100 + 


Figure 1.2: Channel Shifting Approximated by a Smooth Curve. [SEK13] 
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Now, let’s consider reactions where two molecules have to collide in order to 
form some product. These reactions are based on the Law of Mass Action which states: 
In a reaction involving the interaction of two types of molecules, the rate of reaction is 
proportional to the concentrations of the two reactants. [SEK13] 

ky 


The reaction A+B P is such a reaction where the rate of the forward 
k-4 


reaction would be k;|A]-[B] (square braces denote concentrations). To obtain this law, 
it is assumed that the molecules are far enough apart to assume they move indepen- 


dently and, thus, their concentrations are low. [SEK13}] 


Let us look at a ”dimerization” reaction in which two A molecules reversibly 


combine to form a complex C’. The kinetic scheme looks like this: 


ky 


The reaction equations are: 


A 
a = —2k, A? +2k_1C and “ = ky A? —k1C. (2.4) 


Initial conditions at t = 0, are A = Ap and C = 0 where Ap represents the 


concentrations of A at t = 0. Combining the differential equations in (2.4) we have 


In this case A+ 2C is equal to a constant and, thus, equal to Ag. Again, we 
consider the molecules to be in free form. Therefore, the molecules are either of form 
A (per unit volume) or form 2C (per unit volume). The recording of complex values of 
C and given parameters k_;,k,, and Ao leads to a great number of possible graphs to 
evaluate. Thus, a reformulation called ”non-dimensionalizing” is required to reduce the 


number of parameters. This will be discussed later. [SEK13] 


The following biochemical model is used as a component of many models in 


molecular and cellular biology. It is central to the study of enzyme-mediated reactions 


12 


in biochemistry. The kinetic scheme of this model is: 


ky 
E+8_"C 4E+P. (25) 
Roi 2 


Here F is the concentration of an enzyme that catalyzes the transformation 
of the substrate, S (concentration), of the reaction into the product, P. This is ac- 
complished by means of an intermediate enzyme substrate complex (concentration C) 
where the enzyme and complex are bound together. Based on the Law of Mass Action 


and pervious examples, the differential equations for this kinetic scheme are: 


“ = —-k BES +k 1~C+ koC, (2.6a) 
dS 

= E = 2. 
a k, ES + k_1C, (2.6b) 
“ = ki ES = kiC = koC, (2.6c) 
dP 
—_— = ; 2: 
a = kee (2.6d) 


The initial conditions at t = 0 are: E(0) = Epo,S(0) = So,C(0) = Co = 
0, P(0) = Py = 0. There are two conservation statements produced from these differen- 


tial equations and initial conditions. They are: 
E(t)+C(t)=Eo and S(t)+C(t) + P(t) = So. 


These statements indicate that, in both free and bound forms, the total amount 
of enzyme, and the total amount of reactant in the substrate, complex, and product 


forms is constant. Substituting # = Ep — C into the differential equations gives us two 


differential equations with two unknowns: 


< = ki (Eo C)S + k_1C, and 
dC 
Te = ki (Eo C)S (k_1 t k2)C. 


It is often the case that the enzyme-substrate complexes form rapidly, and, 
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thus, the number of complexes are roughly constant. Therefore, a x0. This assump- 


tion is called a quasi-steady state approximation. We can then solve for C in terms of 


S' and Eo to obtain: 


ae where ky», = ee 


hm + 8? ky 
Substituting C and ky, into the differential equation of the product we get: 


dP Vinaze® 
dt = km+S 


where Vinaz = keEo. (2.7) 


This equation is defined as the reaction velocity approximating the rate at 
which substrate is used up and product is formed. It is known as Michaelis-Menten 


kinetics. [SEK13] 


1.2.4 Polymerization Reactions 


Now, let us consider polymerization reactions, in which identical subunits 
(monomers) form polymers that can grow in size. The Law of Mass Action will still 
apply for these reactions. The first polymerization reaction we will look at involves 
simple aggregation. When monomers grow anywhere on the polymer it is described as 
a simple aggregation polymerization reaction. The following Figure 1.3, shows a simple 
aggregation with ky the rate of binding and 6 the rate of disassembly and turnover of 
the polymer. 
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Figure 1.3: Simple Aggregation.[SEK13] 


We define the following variables: 


c(t) denotes the number of monomer subunits in the volume at time t. 
F(t) denotes the amount of polymer (in number of monomer equivalents) at time t. 


( 
A(t) denotes the total amount of material (in number of monomer equivalents) at time 


+ 


We assume the rate of growth is a product of c and F, with rate constant 
ky > 0 and the rate of disassembly or turnover is linearly proportional to the amount of 
polymer, with rate constant 6. The differential equations are then consistent with the 


other models we have looked at thus far and are: 


= = —kpeF +06F, (2.8a) 


dF 
— =krcF — oF. 2.80 

Notice the first terms in these equations describe the association of a monomer 
and a polymer based on mass action. The last terms are the polymer turnover at rate 


6. The total amount A is conserved and for physical relevance, we must have c < A and 
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F <A. Substituting F = A -—c into (2.8a) we get: 


de 6 
—=k;(A-—c)(——c 
r( Ge ) 


dt 


where ~ is defined as the critical concentration of monomers, Cc,it. The right-hand side 


f 
of this equation forms a quadratic equation and its graph is shown in Figure 1.4 below. 


dc/at 


Figure 1.4: Quadratic Equation dC/dt. [SEK13}] 


The arrows on this plot indicate values of c for which c would increase (<2 > 0: 
arrows point to the right) versus places where c would decrease (2 < 0: arrows point 
to the left). We, also, observe stagnate points at which there is no flow. These stagnate 


points are steady states. The case A < C¢piz is shown in Figure 1.5 (below). 


Figure 1.5: Steady State of Equation dc/dt. [SEK13] 


16 


Since c < A, part of the state space is blocked and is represented by the “gray 
zones”. The diagram summarizes the qualitative behavior of the model. Notice, also, 
that a nontrivial level of polymer occurs only if A > Cc;jz and that steady state occurs 

dc é 


when 4 = 0. Thus, the amount of monomer left in solution is either Cori, = kp OF A. 


Therefore, the amount of polymer, F’, at steady state is equal to zero. [SEK13] 


Now, we will look at polymers with growth only at their tips (end of their 


filaments) as shown in Figure 1.6. 


monomer 


Figure 1.6: Polymer Growth at Their Tips. [SEK13] 


In this case, we define n as the number of filament tips at which polymerization 
can occur. We consider this number to be constant, and that breakage and branching 


do not occur. Thus, the model is: 


“ =-kren+o6F, (2.9a) 
F 
“ =ksen—d6F. (2.90) 


Here, monomer addition occurs at a rate proportional to n. As before, conser- 
vation holds and replacing F' with (A — c) into (2.9a) leads to: 
dc 


a 5A — c(kpn +0). 
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In this case, there is only one steady state with monomer level 


5A 5 
=~" _=6A where B= —°—. 
a ee a 


In addition, the factor 6 satisfies 8 < 1 since n, kr > 0. That means that this 
steady state exists in all cases, unlike the previous situation where up to two steady 
states were possible. Since the number of tips, n, is constant, the steady state levels of 
monomer and polymer are proportional to each other, so that adding monomer to the 
solution will increase both forms. We also see that the growth of the polymer is linear. 
Thus af = constant. Figure 1.7 distinguishes the differences of (a) the case of a simple 
aggregation for monomer concentration Co,;, <A, where no polymerization occurs and 


(b) the case of growth only at filament tips. [SEK13] 


F. 


(a) 


(b) 


Figure 1.7: Differences of Polymerization Growth. [SEK13] 
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All the modeling techniques discussed above have illustrated ways of describ- 
ing rates of simple chemical reactions. These modeling techniques provide a means for 
reaching reasonable approximations and for providing estimates for behavior under var- 
ious situations. At this point, we can see and appreciate the usefulness of mathematical 


modeling in a biological environment. [SEK13}] 


1.3. Non-Dimensionalization and Scaling 


Now, we will look at non-dimensionalization and scaling. We can check the 
consistency of model equations by reformulating a model in terms of dimensionless quan- 
tities. This reformulation will ensure that all terms have the same set of units in an 
algebraic or differential equation. Non-dimensionalizing a model reduces the number of 
free parameters and reveals a smaller set of qualities that govern the dynamics. Deter- 
mining larger or smaller magnitudes helps approximate solutions when using asymptotic 
analysis techniques. We will show that dimensionless variables can be selected in a va- 
riety of ways using the concept of scaling. We will use several of the kinetics problems 


previously discussed to show this procedure. [SEK13] 


1.3.1 Kinetic Scheme 


ky 
Let us consider the kinetic scheme A B and the derived equation 
ky 
dA 
Th = —(kj +k_1)A+k_1M, A(0)=Ap and M=Ap+ Bo. (2.10) 


Remembering rate constants k; and k_1 have dimensions of 1/time, both dA 
and (k; +k_,)A must have the same dimensions, concentration/time. For definiteness, 


let us define a dimensionless time, t*, by 


t 


ety 


ie. th =k yt. (2.11) 


Note: *’s will denote variables carrying no dimensions. The natural way to non- 
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dimensionalize A is via its initial concentration. Thus, we define a dimensionalized 
concentration 

A 
—., 2.12 

= (2) 
Using the chain rule and adopting (2.11) and (2.12) we obtain 


dA d(Apa*) da® da* dt* da* 
as Se ee ay ea eee 
dt dt Ov dt Catt dt Oat 


k_1. (2.13a) 


or by direct substitution 


dA d(Aga*) da* 
dt d(t*/k_-4) Bat0 gg SE) 
Therefore, 
da* ° 
kiAoT = —((ky + k_1)Ao)a +k 1M. 


Now dividing both sides of the equation by Agk_, we arrive at 


da* Ao k_41M 
= L fo_ ‘i ‘ 
(Ai 1 Aga ae 


da* _ Ae) eA M pu M 
= ea ee 


We can define two dimensionless parameters, 


ga M _ Act Bo 
ky” ~ Ap Ay 


Then the resulting dimensionless equation is 


da* 1 
= —-a* —a*. 2.14 
or a0 +0-a ( ) 


Furthermore, the initial condition A(0) = Ap leads to a*(0) = A(0)/Ap = 1. 


20 
Thus dropping the *’s, the new version of the model is 


= a+@—-—a, a(0)=1. 
dt a* yO) 
where € and @ are dimensionless parameters. Clearly, non-dimensionalizing and scaling 


the original problem decreased the four dimensional parameter model to a two non- 
dimensional parameter model. [SEK13] 


1.3.2 Dimerization Model 


Now, recall the dimerization model, 


A 
a = —2k,A*4+2k_1C, (2.15a) 
“ = k, A? — k_1C. (2.15b) 


ky 
at t= 0, A = Ap, and C = 0 with kinetic scheme A+ A C 
eal 


Similar to the previous example, we define the dimensionless time and concen- 
tration variables as follows: 


t kat » A x C€ 
ee ae Sa CSS 
ey aa Ao 
Therefore, 
Le ¥ * 
t= ; A= Aga’, C = Aoc’, 
k_4 


(2.16) 


and the initial conditions can be written as t* = 0, 


a*=1, and c*=0. 
Substituting (2.16) into (2.15a) and simplifying, we get 


dA d(Aoa*) 


- = —2k;(Apa*)? + 2k_1(Aoc*). 
dt d(t*/k-1) AS ea ae) 


This leads to 
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da 


_1A 
ky 0 


= —2k A2(a*)? + 2k_1Apc*. 
Then simplifying, we get 


da* ky Ao 
dt* k_4 0 


Similarly, substitution of (2.16) into (2.15b) yields 


= —20(a*)? +2c* with 0= 


dc* 
dt* 


=a") = 


Dropping the *’s, the new dimerization model becomes 


“ =-20a7+2c, a(0)=1 (2.17a) 
d 
FT =6a?—c, (0) =0. (2.170) 
Note 6 = Miao is dimensionless since the concentrations, denoted [ ], of Ao, ki, and k_; 


are [Ao] = L~3, [k_1] = T~1, and [k,] = L°T~!. The dimensions of [k;] are L°T~! be- 
cause the dimensions of k; A? and k_1C in equation (2.15) must be the same. Therefore, 


[k,] = L3T—! and @ is dimensionless.[SEK13] 


We observe from the new model that the dimensionless complex concentration 
c depends only on the single dimensionless parameter 9. Therefore, the variable com- 
bination C/Ao is not a general function of the parameters k_1,k,, Ao and time t but 
rather a function only of the combinations (kj Ao)/k_1 and k_,t. Thus, all the data 
concerning dimerization can, in principle, be presented on a single graph of c = C/Ao 
as a function of t* = k_,t for a number of different values of 0 = (ki Ao)/k_1. See Figure 


1.8 below. [SEK13] 


22 


10 a 
rl 
lA _—_——— $=P 
5 / Pe 
| / po o=¢1 
- 
f 1 = —> kif 
0 l 2 3 


Figure 1.8: Dimerization Data (C = C/A,). [SEK13] 


1.3.3. Polymerization by Monomer Aggregation 


Lastly, let us create a dimensionless model for polymerization by monomer 


aggregation. Recall equation (2.8c) in the form 


= (—kyc+8)F =(-kyc+6)(A—c). (2.18) 
Rescale time by (1/6) and rescale concentration by the total amount A (both 
of which are constant). This means substitute t = t*/d and c = c*A into equation 
(2.18), then using substitution as before, 
d(c* A) dc* 


Ge]e ~ Oo age = (hrcA + 8A cA), 


Simplifying and multiplying both sides by the constant (1/A6) leads to 


de® 1 eo #4) (IA : 
ae Tag hset A + 6)(A— eA) = (- et +:1)(1 — e*). 


Dropping the *’s we arrive at 


dc A 
=(l-—ac)(l1—c) where n= eg 


which is the ratio of total amount A to critical concentration Cpt = 6/ ky and dimen- 


sionless. Thus, the inherent structure of the model depends only on the grouping of 


parameters in a. [SEK13] 
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The new equations developed in these examples are written in their simplest 
form, containing only dimensionless variables and parameters, and were obtained with- 
out solving any equations. The number of dimensionless parameters is generally smaller 
than the original number of parameters which makes theoretical manipulations easier. 
Thus, conclusions developed from the reduction in parameters can be of great impor- 
tance for experimental or numerical work. This minimizes the amount of experimen- 
tation that is necessary to describe the possible variation of the results for different 


parameters values. [SEK13] 


1.4 Geometric and Qualitative Methods to First-order Dif- 


ferential Equations 


Now, we will apply qualitative and geometric methods to first-order differential 
equations, where the goal is to construct sketches of the solutions that do not require as 
much technical work. We will introduce the role of parameter sensitivity and illustrate 


some transitions, called bifurcations, that take place as a parameter is varied. 


1.4.1 Stability of Steady States 


To begin, let us first understand the stability of steady states of a first-order 


differential equation. Consider the general differential equation 


dx 

= 2.19 
T= f(a), (219) 
with a steady state value x = x5. It follows that 


akss 
dt 


=O. fies) = 0: (2.20) 


Let x(t) = ass + 2p, where x» = xp(t) is a small-time dependent quantity, 
called a perturbation, of the steady state. When z,(t) decreases with time, we see that 


x(t) returns to the steady state value and we say that the steady state, x55, is locally 
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stable. When x,(t) increases with time, we see that «(t) moves away from the steady 


state value and we say the steady state, xss, is locally unstable. [SEK13] 


To see what truly happens, let us substitute x(t) = 25,;+2, into the differential 


equation (2.19) as follows, 


dla ss a Lp] 


dt om igs + Le) 


When we simplify the LHS using the additive property of derivatives and the 


RHS by using the Taylor series approximation since 2, is small, we get 


d|%ss+2p] diss Zz dip 
dt dt dt 


2 
L 4b; " 
~ f (Zs) a5 Lyf (ages) + oe (ass) + h.o.t., 


where h.o.t stands for “higher-order terms”. Using the steady state equations of (2.20), 
we can eliminate some of the terms from each side of the equation to obtain 

dx , x2 ” 

—? & apf (rss) + a (%ss) + h.o.t.. 

In addition, the magnitude of terms that are quadratic or of higher power in 
small perturbations are very small since, as the powers increase, the magnitude of the 
quantities decrease. Thus, cs = 0 and it = 0. Therefore, we obtain the following linear 
equation governing the perturbations: 

dip 


WE wm f (Lss)Xp, 


where f (x55) is a constant of negative, positive, or zero value. Denoting f (x55) by A, 


we have 


which has the general solution 


x(t) = Cexp(At). 


Thus, the linear stability condition is: 
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=f’ (ass) > 0 => exponentially growing perturbations, then x,, is unstable. 
A= ‘i (ass) < 0 = exponentially decaying perturbations, then x5 is stable. 
A=f' (x55) = 0 => no conclusion, since higher-order terms are needed to determine 


local behavior near the steady state. [SEK13] 


For reasons of analogy to higher-order systems of differential equations, the 
quantity 4 = f' (ass) is called an eigenvalue of the ODE at the given steady state. It 
is seen that there is only one eigenvalue, whose value is a real number, at any steady 
state of a single first-order ODE. In such cases, the eigenvalue can be interpreted as a 
rate of growth (if positive) or rate of decay (if negative) of small deviations from the 
steady state. We can restate the stability results verbally by saying that for a single 
ODE, stability (or instability) of a steady state is equivalent to finding the eigenvalue 
is negative (positive). A zero eigenvalue is neutral and is usually a sign that some tran- 


sition in stability is at hand. [SEK13] 


Now, let us consider an example where there are three steady states 


d. 1 
x 53 


A C(x 3 )= f(x), c>0 constant. (2.21) 


Solving for the steady states (2 = f(x) = 0), we find that there are three such points, 
one at x = 0 and the others at x = +3 and x = —\V/3. These are the intersections of 


the cubic curve with the z-axis in Figure 1.9a below. 
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Figure 1.9: (a) Steady-State Points. (b) Steady-State Flow. [SEK13] 


We surmise the direction of flow from the sign of f(x), and use that sketch 
to conclude that 2 = 0 is unstable while both x = —V3 and x = +V3 are stable. In 
Figure 1.9b, we show numerically computed solutions to this equation with a variety 
of initial conditions. We see that all positive initial conditions converge to the steady 
state at « = +/3 = 1.73, whereas those with negative initial conditions converge to 
x = —V/3 & —1.73. Thus, as seen, the outcome depends on the initial conditions and is 
an example of bistable kinetics. Figure 1.10 is an abbreviated way of representing the 
solutions to our equation and is a sketch of the flow along the x-axis called a “phase 


line”. [SEK13] 


wl 


Figure 1.10: Phase-Line. [SEK13] 


1.4.2 Bifurcations 


Now, let us look at the same differential equation, but introduce a new param- 
eter and explore some of the behavioral transitions (bifurcations) as that parameter is 


varied. 


“ =o(n— 32° +A)= f(a), (2.22) 
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where A is some additive positive or negative constant. We can set c = 1, without loss 


of generality since time can be rescaled. 


Now, look at Figure 1.11la below. 


ix) 
Nx) f(x) 
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Figure 1.11: (a) Curve Shift. (b) Flow of Shift Curve. [SEK13] 


As we know, A shifts the curve up when A > 0 and down when A < 0. As 
seen on Figure 1.11b, when A changes, so do the positions and number of intersection 
points of the cubic and x-axis. The black dots indicate stability while the white dots 
indicate instability. [SEK13] 


There are large values of A in both the positive and negative directions beyond 
which two steady states coalesce and disappear. At each of these values the graph of 
f(x) is tangent to the z-axis and is shown below in Figure 1.12. This type of change 
in qualitative behavior is called a bifurcation and A is called a bifurcation parameter. 
[SEK13] 
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Figure 1.12: Change in Quality Behavior of Bifurcation Parameter. [SEK13] 


The behavior of an entire system can be summarized by assembling a bifur- 
cation diagram with the number and relative position of its steady states (or more 
complicated attractors) along the y-axis and the variation of the bifurcation parameter 
along the z-axis. Figure 1.13 represents the bifurcation diagram of equation (2.20). 
[SEK13] 


Figure 1.13: Bifurcation of Equation 2.20. [SEK13] 


Figure 1.13a simply rotates Figure 1.11b and removes the arrows giving a clear 
view of the steady states along the y-axis while Figure 1.13b is the true bifurcation di- 
agram where the dotted line represents the unstable state (white dots). Because the 


bifurcation curve appears to fold over itself, it is known as a fold bifurcation and has 
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two bifurcation points: one at a positive point and one at a negative point. [SEK13] 


The existence of two stable states in a differential equation model is often 
described by the term bistability. This behavior occurs in many biological situations. 


Bistability is accompanied by the following hysteresis. [SEK13] 


x 


Figure 1.14: Hysteresis of Two Stable States. [SEK13] 


Let the parameter A start as the negative steady state value. As A gradually 
increases we remain at steady state, but the value of that steady state disappears and a 
rapid transition to the positive steady state value takes place. Now, if we let the value 
A decrease back to lower values, the elevated state moves left along the upper branch 
until the negative bifurcation value of A is reached. This type of hysteresis is often 
used as an experimental hallmark of multiple stable states and bistability in a biological 


system. [SEK13] 


To generalize, let us consider the single differential equation 


dx 


dt at fr(z), 


where f depends on one parameter “r’ that takes the role of the “bifurcation param- 
eter”. First consider the condition f,(x) = 0. This simply restricts our attention to 


steady states of the ODE. 
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As r varies, the relation f;(2) = 0 corresponds to a set of curves in the ra 
plane (the bifurcation plot) that are smooth functions of r except at special points 
where f.(z) = 0. A mathematical result states that these special points are the only 
places where branches of steady states can come together (bifurcate). A point in the rx 
plane satisfying both f,(x) = 0 and f.(2) = 0 is a bifurcation point, and the value of 
the parameter r at such a point is the bifurcation value. The conditions for bifurcation 
are hence 


fr(x) = 0, f(x) =0, Bie 2S eae Pf SO: 


Geometrically, these conditions indicate that the function f,(x) intersects the 
x-axis and is tangent at that point. Analytically, these conditions imply that at the 


steady state x = rs, there is a zero eigenvalue. [SEK13] 


1.5 Mathematical Model of a Notable Biological Problem 


Now, we will use the tools we have developed to construct, analyze, and inter- 
pret a mathematical model of a notable biological problem. As a case study, we take 
the spread of an infection in a population of initially healthy individuals. Using di- 
mensional analysis, qualitative methods, and ideas of bifurcations, we will consider how 


parameters affect the behavior of the model and what this implies biologically. [SEK13] 


1.5.1 Derivation of Model 


First, we will derive a model for the spread of infection. Assume any individual 
is equally likely to come into contact with any other individual and subdivide the popu- 
lation into two classes: those that are healthy, but susceptible to the disease, and those 
that are infected and able to spread the disease through contact. Figure 1.15 illustrates 
a view of the process that is closed (the population neither increases nor decreases). 


[SEK13] 
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Figure 1.15: “Closed” Process for Speed of Infection. [SEK13] 


The “reaction scheme” for this illustrated view is: 


S47-397, J—8. (9.93) 


To simplify, we will also assume that infected individuals recover at a constant 
rate and become susceptible again with no immune period. Let us define S(t) = the 
number of susceptible people and I(t) = the number of infected people in the popula- 
tion. The spread of infection requires contact between healthy and sick individuals, as 
indicated in the reaction scheme (2.23), the rate at which this type of contact occurs 
can be approximately represented by the Law of Mass Action. Thus, the rate of in- 
crease of infected individuals would be proportional to the product SJ. Let us call the 
proportionality constant 6. Since the rate of recovery of a sick individual is assumed to 
be constant, the overall rate of “flow” out of class J and into class S is proportional to 


I. Let us denote the rate by uw. [SEK13] Then the equations are: 


dS 
aE pl — BSI, (2.24a) 
c = BSI — pl. (2.24b) 


Note the flow into one class is identical with flow out of the other class. Therefore, the 


total population, N(t) = S(t) + I(t), is conserved as shown: 


dN dS di 
ca a ee 


Therefore, NV is a constant. By conservation, we can eliminate J(t) from the equations 


and substitute I(t) = N—S(t). Since N is constant, if we find S(t), then this relationship 
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provides I(t). Thus, it suffices to keep only one equation. We will keep equation (2.24a), 


rewritten as 
as _ 
dt 
together with I(t) = N — S(t). The model is thus reduced to a single ODE in S(t). 


u(N — 8)— BS(N —S), (2.25) 


Equivalently, we could choose to eliminate S rather than J. This would lead to 


dI 


OT O(N - DI al. (2.26) 


1.5.2 Dimensional Analysis and Scaling 


Next, we will apply dimensional analysis to the model to reduce the number 
of parameters and make the model dimensionless. Let us measure time in units of days, 
then the left-hand sides of Equations (2.24a) and (2.24b) have dimensions of [number 
of people] / [time]. Therefore, ~ must have units of 1/time and @ should have units 
of per person per unit time. Since yz has units of 1/time, it follows that 1/p carries 
units of time. This time unit is a typical time associated with recovery from infection. 
Let «*(t) be defined as the fraction of the population in the infected class and y*(t) be 
defined as the fraction of the total population in the susceptible class. [SEK13] 


It is convenient to define dimensionless variables: 


oT] N’ v N’ =a 
Since S(t) + I(t) = N = constant, it follows that 
ee ee ee Saas 
Ny ee 


Rewriting the relationships obtained in a form that can be used to substitute 


directly into the model equations of (2.24), namely, S = y*N, I = a*N, t = t*/wp yields 


= px" N — B(y* N)(2*N), (2.27a) 


——_~ = B(y*N)(2*N) — pa*N. (2.27b) 
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Canceling the constant common factors of N and yu from both sides we arrive 


at 
ae eGR Nes 
ae ( 5% ara (2.28a) 
da* BN 
= *y ae 222 
east -2 (2.280) 


Note (&) is the single remaining ratio of parameters that we will denote as 
Ro. Ro is an important quantity since its parameter combination governs qualitative 


behavior. Rewriting the equations in terms of R, and dropping the stars leads to 


4 = — Ro xy; (2.29a) 
dx 

—=Rh, — 2. 2.2 

a Ro xry-—2x (2.295) 


Let us illustrate the process of scaling from the perspective of the reduced 
model (2.26) where conservation, y* + 2* = 1, was used. Substitutions of J = «*N and 
t = t*/p into this model leads to 

d(a*N) 


eye) = B(N —a*N)a*N — pa*N. (2.30) 


Canceling factors N and yu from both sides and dropping the stars leads to 


\1-a)r-2, (2.31) 


which together with y = 1—xz, completely specifies the problem. The same dimensionless 


parameter ratio Ry = oN is seen here. We can rewrite (2.31) as 


dx 
—=R,(1-a2)x—-«z. 
a = Ro(1 —2) 

Considering the fact that 1/y is a typical recovery time, when there is a single 
infected person, J = 1, and the population of susceptibles is S ~ N, then the number 
of new infections per unit time is 8N x 1 = BN. Now, 1/, is the typical time that the 
infected person is ill and can spread the disease. Thus, the total number of new infec- 


tions stemming from this single infected individual is oe = Ro. The development of 
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the model is complete. Our next task is to analyze and understand its behavior. [SEK13] 


Let us write a single ODE for the infected fraction of the population in several 
suggestive forms, namely 
dx 
ae R(l-s)r-xc=2(Roy-1)=2((Ro-1)- Roz]. (2.32) 
1.5.3. Steady State 


Then the steady states of the model are 
0=2[(R, — 1) — Roa], (2.33) 


either « = 0 and then (by conservation) y = 1. This corresponds to a population 
that has no infected individuals, I(t) = 0, S(t) = N. We will refer to this as 
the disease-free equilibrium. A second possibility is that z[(R, — 1) — Rox] = 0 so 
x =(R,—-1)/Ro =1-—(1/R.). Using conservation once more, we find that in this case 
y = 1/R,. This steady state has some proportion of the population in the infected class, 
and is denoted the disease endemic state. However, we observe that this steady state 
is biologically feasible only if x > 0, which means that R, > 1. When this steady state 
exists (which implies that it is positive), we say that the disease can become endemic, 


which means that it can take hold of some constant fraction of the population. [SEK13] 


To summarize, 


1 1 
Disease free: x1, = 0, Yo =1, Disease endemic: x; = 1 n= 2 (234) 
Ro Ro 


To convert these findings to results for the unit-carrying variables, we multiply 


each quantity by N, so that S = yN, I = <zJN, and thus 


Disease free: I, = Nx, =0, So=N, yo = N, (2.35a) 


1 1 
Disease endemic: I; = Na, = N(1— Rp” 5, =Ny = Ne: (2.35) 


‘O O 


This result is summarized in the following Threshold Theorem: 
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Theorem 2. In a simple SJ disease-dynamics model, the disease can become endemic 


BN 


only if Ro > 1, where R, is the reproductive number of the disease, and Rp = i 


[SEK13] 


1.5.4 Qualitative Techniques - Behavior of the Model 


We now apply qualitative techniques to understanding the behavior of the 
model. We will use the ODE for the infected fraction of the population (2.32) but 


written in a more convenient form 


dz 1 
— = Ror/1— —)—2]= 2.36 
7a Botlt-Z)-a1=f@) 2.36) 
Note that the expression (1— ze) is one of the steady state values obtained in 
(2.34), a constant. Figure 1.16 shows the flow diagram for the RHS of f(x) in (2.36). 


Se) (a) Je) b) 


1- IR 


4 


0 x x 0 


Figure 1.16: Flow Diagram for Equation 2.36. [SEK13] 


We see that in the case R, > 1, the disease will progress towards the endemic 
steady state, x = 1— (a and in the case R, < 1, the only steady state is at x, = 0, so 
the disease is eradicated. Stable steady states (black dots) have arrows directed towards 
them, and unstable steady states (open dots) have arrows directed away from them. The 
simulation of this model is shown below in Figure 1.15. It confirms our results that all 


initial conditions with x > 0 eventually approach the endemic steady state. 
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Figure 1.17: Simulation of Model for Equation 2.36. [SEK13] 


Now, we will confirm the linear stability of those steady states. Let us recall 
that for an equation of the form 4 = f(x), a steady state x,5 (satisfying f(x) = 0) is 
stable whenever f (x55) <0, and unstable when f'(a,5) > 0. 


For the problem at hand, we have 


Hence, 


f (2) = Rol( - 3) ~ 2a 


Thus, for the disease-free steady state, r, = 0, we have 
f (x) =f (0) = R,[(1 - —)] = Ro -1>0 when R, > 1. 


This means that the disease-free state is unstable when R, > 1 (and stable 


otherwise). Similarly, for the disease-endemic steady state, 7,, = 71 = 1 — re so 


1 1 1 
I] = -Rol1 — 


f (x1) = Rol(1 a 


Thus, the disease-endemic steady state is stable only when R, > 1. [SEK13] 


We now develop an explicit relationship between the size of the parameter, Ro, 


and the steady state level in a one-parameter bifurcation diagram. 
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We have two steady states. One of these, 75 = 0, does not depend on Ry. The 
second appears when R, = 1 and is then 2; = 1— (a As R, — ov, this second steady 


state will approach 1 as shown in the following Figure 1.18. 
x 


06 


04 


Figure 1.18: Bifurcation of Model for Equation 2.36. [SEK13] 


The solid line corresponds to the stable steady state and the dotted line to the 
unstable steady state. This type of bifurcation is a transcritical bifurcation, in that, 
for most values of R, there are two steady states, but for R, = 1 these steady states 
merge and exchange stability, whereas with the “fold” bifurcation they merged and dis- 


appeared. [SEK13] 


The analysis, simulation, and bifurcation plot show that the outcome of the 
infection depends on a single parameter Ry, = a and that there is a transition in the 
qualitative behavior at Rj, = 1. Below R, = 1, the disease cannot “reproduce” itself 
fast enough to be sustained and consequently disappears after some transient. Above 
R, = 1, this “reproductive number” implies that each infected person infects more than 


one uninfected person, on average, and the disease becomes endemic. [SEK13] 
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Chapter 2 


Practical Applications 


The following applications are to test our understanding of everything we have 
learned thus far. We have taken these applications directly from the book, ” A Primer 
on Mathematical Models in Biology” by Lee A. Segel and Leah Edelstein-Keshet pages 
111-114. 


2.1 Application 1 - Total Population 


a) Show that eliminating the susceptible variable from the system of equations gives 
rise to equation 2.26. 

b) Simplify that equation and find its dimensionless version using the same procedure 
as shown in the examples. 

ANSWER: 

a) Total population is N(t) = S(t) — I(t). Therefore, S(t) = N — I(t) (N is constant) 
and substituting this equation into a = BSI — pl gives a = B(N—-—I)I- pl. 


b) Time will be measured in units of days with a and a having dimensions of the 


A 
tame 


number of people/time, then jz has units of and ( has units per person per unit 
time. Defining x*(t) as the fraction of the population in the infected class and defining 
y*(t) as the fraction of the population in the susceptible class, we define dimensionless 


variables 


Since S(t) + I(t) = N = constant, it follows 


ee ee Sea 
Ma oe 
Now, substituting, S = y*N,J =2*N,t= 7 into (2.24a) and (2.24b) we get 

d(y* N) 

OY ON! = ye*N — B(y*N)(2*N), 

Fegy 7 BEN Btu VED) 

d(a*N) 

——_~ = B(y* N)(a*N) — pa* N. 

dey, ee ee 


Dividing both sides by N and yu leads to 


dt* Ls 


BN 


Rewriting the equations with R, = Carve and dropping the stars, we get 


= =«— Rory, 


« = Rory — &. 


Using y* + 2* = 1 and substituting J = «* N and t = t*/p into 


dI 


—=P(N—-I)I-pl 
rats )\I- pl, 
gives 
d(a*N) 
= B(N —a*N)a*N — pa* N, 
ae]n) aie 
leading to 
dx BN 
— 1 
Ba )a- 20-2 


Then, we divide by and N with R, = (=) to get 


39 


40 


dx 
ae Ro(1-—2)x — «2. 


Using y* + 2* = 1 and substituting S = y*N and t = t*/p into 


dS 
Ge MIN — S) — BSN — $), 
gives 
d(y* N 
ean = UN —y*N)— B(y*N)(N — y*N), 
leading to 
dy _ BN 
oe 7 y1—y) 


Then, we divide by and N with R, = (=) to get 


“ = (1—y) — Roy(1 — y). 


2.2 Application 2 - Find Steady-State 


Find steady states corresponding to the endemic disease in the SJ model in 
terms of the original, dimension - carrying variables (rather than the dimensionless vari- 


ables x, y). [SEK13] 
ANSWER: 


We do not have to redo the work - merely “convert” back from dimensionless 


to unit-carrying variables. The disease endemic non-dimensioned steady state is: 


To convert to dimension-carrying variables, we multiply each quantity by N: 


N N 
Nee. INE 
ot Fe NE pe 


substituting 6N/p for Ro 

LN 
— —— Ny = 5S, = >=. 
BN? ES PL BN 


Thus, the steady state of the original, dimension-carrying variables is 


2.3. Application 3 - Endemic Intervention 


Al 


Suppose an emergent disease threatens to become endemic. Based on the anal- 


ysis in this chapter, explain what interventions might be used to suppress it. Use Ro 


or parameters of the model to explain how various interventions affect the dynamics. 


[SEK13] 


ANSWER: 


N 


Ro = BN, where GN is the number of new infections per unit time and 7 is 


the typical time that the infected person is ill and can spread the disease. 


From the stability analysis, we know that the steady states are at Ro = 0 and 


R, = 1. If we maintain the value of R, between 0 and 1, then the disease cannot sustain 


itself and, thus, will not spread. To keep R, between 0 and 1, the value of 6N must 


remain less than y. Therefore, we want 1/j to decrease by using a variety of medical 


care or a pattern of medical care that can decrease the time an infected person is ill and 


can spread the disease. 


Also, we want the number of new infections per unit time to decrease. Thus, 


isolating the infected patients, decreasing the amount of time an infected person is in 


contact with an uninfected person (by awareness of symptoms and fast action), and 


communicating with the public to learn the facts of the disease and avoid infection are 


critical. 
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2.4 Application 4 - Disease Outbreak Prediction 


A United Nations report about a recent disaster in a third-world country pre- 
dicted that as refugees crowd into relief centers, disease outbreaks might occur. Does 


the model support this assertion? Why or why not? [SEK13] 


ANSWER: 


Yes, the model does support this assertion. As discussed in Application 3, Ry = eal 


pw? 
and, therefore, to keep disease from becoming endemic @N must remain small. Since 
BN is the number of new infections per unit time, increasing the potential of exposure 
by bringing a large number of people to relief centers increases the rate of possible 


infection by mass action. Since the rate of mass action is proportional to the product 


STI and,thus, 3, this would increase BN. 


2.5 Application 5 - Vaccination 


Vaccination is sometimes effective at preventing epidemics. Here we will sup- 
pose that vaccinating a fraction p of the population is equivalent to protecting pN 
people from being either susceptible or infected. This effectively reduces the size of the 
population that can ”participate” in the disease dynamics. Determine the fraction p 
that would have to be vaccinated in each case to prevent an endemic disease. [SEK13] 
a) Smallpox, for which R, ¥ 5. 

b) Polio, for which R, ¥ 6. 

c) Measles, for which Ro © 12. 

ANSWER: 

a) ke = oN is the total number of new infections from a single infected person. The 


steady state Disease Endemic is 


r=l 


with x representing the fraction of the population in the infected class. Thus, with 


Ro & 5 we get 
1 : TL 2.8 : 2 
1 5 ’ Y1 5 
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Therefore, to prevent a smallpox endemic, we need to prevent the population 
from getting to the steady state of 7, = .8, so the fraction of the population that needs 


to be vaccinated is p = .2 or at least 20 percent. 


b) With Ro ¥ 6 the steady state Disease Endemic is 
1 1 
gy=1- eo 1 — .16667 = .83334, y= 5 .16667. 


Therefore, to prevent a Polio endemic, we need to prevent the population from getting 
to the steady state of 7, = .83334, so the fraction of the population that needs to be 
vaccinated is p = .16667 or at least 16.67 percent. 

c) With R, © 12 the steady state Disease Endemic is 


1 1 
Ly =-1— Tia 1 — .08334 = .91667, y= aa 08334. 


Therefore, to prevent a Measles endemic, we need to prevent the population 
from getting to the steady state of x; = .91667, so the fraction of the population that 
needs to be vaccinated is p = .08334 or at least 8.34 percent. 


2.6 Application 6 - Separation of Variables 


Use separation of variables to find an analytical solution to Equation 2.36 with 


initial condition x(0) = a, where 0 < « < 1— os It is advisable to first recast the 


equation by defining the constant 


to obtain 
o = Rea[f 2], O< 2 <.B: 


Compare your solution to the results described in this chapter. [SEK13] 
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ANSWER: Equation 2.36 is 


dx al 
= Retl(1- =) - 2] =f) 
Substituting 6 = (1 — re we get 
dx 
= Ro«(6 — 2) 


To integrate, we will use partial fraction decomposition. Thus, we have the partial 


fraction equation as 
1 A B 


Rox(B — 2) ~ Rox 
Dividing both sides by Rox(3 — x) leads to 


l= A(6 _ L) + BRox. 


Letting « = 0, we can solve for A as follows 


1= Ag, 
1 
A= =, 
B 
Therefore, our equation becomes 
—2£ 


Letting « = 1, we can solve for B as follows 


—1 
1 = “3 + BRo, 
B-1 B-B+1 
Ro Ro BRo 
Therefore, our original equation becomes 
1 1 1 


Roa(B—«)  BRot  BRo(B—2) 


45 


Our integration will then be 


| awacy”- | ee” + | eae - | 


which is 


Ina In(6 — x) +C =t, where C is a constant. 


Since x and (3 — x) represent numbers of people, they are positive and can be 


written without absolute value bars. Multiplying this equation by GR, we get 


Inz —In(8 — x) + BRoC = PRot. 


Now, using the properties of logarithms and raising each term to the value of 


e, we get 
oC 

gor cbRot 

(8-2 
Solving for x gives us 

BePRot 
C= x(t) => GBRoC 4. cBRot a cBRot’ 
with x(0) = 2, we can find x, as 
pee B 


10.53 GBReOrt 60 =, gona | 


Thus, 
ePRoC — Bete 
Lo 


p] 


and finally, we obtain the analytical solution to equation 2.36 as 


BaoeP Ret 
B= Bet tgeh hor 


x(t) = 


2.7 Application 7 - Susceptible/Infected Pool Model 


Suppose that the population has births at rate b to the susceptible pool and 
mortality at rate 6 from the infected pool as shown in the schematic diagram of Figure 


2.1. [SEK13] 
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a) Write down the modified equations of the model. 

b) Determine if conservation holds in this case. 

c) Determine a dimensionless form of the model using the same kind of method as em- 
ployed for the example in this chapter. 

d) Write an XPP file to simulate this model and plot some of its solutions. 


ANSWER: 


B 
St) Lo] Ie) 
UL 
I ‘ 
Birth Death 


Figure 2.1: Schematic Diagram of the Population with Births and Mortality Rates 


a) Modified equations of the model are: 


dS 
=f = 6ST £55, 
ae BSI +S. 
dl 
— = 6SI-ypl—dl. 
qo Pete 
b) Conservation: Total population is N(t) = S(t) + I(t). 
Thus 
dN dS dl 
= = pl — BSI SI -pl —6Il =bS — 61 
7 Ae a pl — BSI+bS+ 6 [i ; 


which is not equal to zero. Therefore, N is not constant and conservation does not hold. 
c) The dimension of a and we is still [pumber of people] / [time], uz still has units of 
1/time, and (,b,6 have units of per person per time. Defining x*(t) and y*(t) as before, 


the dimensionless variables are 


Thus, 


* 


t 
S=Ny*,1l=Na*,t=-—. 
mM 


Substituting these into a) we get 


Multiplying each term by Na we get 


di ca 5 BN aes, 
dt* = m y T tie ’ 
dt* ps i 


dy — BN b 
dx PN) 5 


d) XPP File, ” disease.ode’: 

S’ =mux*I—betax S*I+b%*S 

I’ = —mu* I +beta* S « I — delta I 

par mu = 0.1, beta = 0.2, b = 0.05, delta = 0.05 

init S=1,1 =0 

@ xp = S, yp = I, xlo = 0, xhi = 2, ylo = —0.2, yhi = 1.2 


done 


A simulation of this model with some of its solutions is: 


AT 
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1 t 
08 + - 
0.6 + \° 

b 
04 F a. 
02 + — 
d 

0 i= 

02 Ls 4 4 4 4 1 4 1 4 
0 02 04 0.6 08 1 12 14 16 18 2 


Figure 2.2: XPP Simulation with Curve a) S=0.2,I=0.8 Curve b) S=0.5,I=0.5 Curve 
c) S=0.8,I=0.2 Curve d) S=1,I=0 


2.8 Application 8 - Kermack and Mckendrick Epidemic 


Consider the simple model for an epidemic due to Kermack and Mckendrick 
[76] discussed in Brauer and Castillo-Chavez [15]: 
dS dl dR 


Go Pele eel a, 


Here R(t) denotes a removed (either immune or dead) class of individuals. [SEK13} 

a) Interpret the equations and sketch a schematic diagram for this model. 

b) Explain why the model can be studied as a 2-variable model; which variable is not 
coupled to the other two. 

c) Consider the “trick” of dividing ue by as Then 


dI/dt _dI _ BSI—4l 
dS/dt dS —BSI ~ 


AQ 


Simplify this equation to obtain an ODE for J as a function of S. Show that you can 


integrate this to obtain the solution curve 


I(S) =—-S 4 3M) + ky, where K is constant. 


ANSWER: 
a) 


s(t) / "| mt) - “~ + Rit) 


Figure 2.3: Schematic Diagram of Simple Model for an Epidemic due to Kermack and 
Mckendrick 


Flows inward to the block contribute positively to the rate of change, whereas flows 
outward contribute negatively. 

The equation dS/dt = —GSI is the rate of change when a susceptible individual 
becomes infected. 

The equation dR/dt = —yI is the rate of change when an infected individual becomes 
immune or dies. 

The equation dI/dt = BSI -— pI is the rate of change between dS/dt and dR/dt. 

b) Since dR/dt is made up of immune or dead individuals, they can no longer be in the 
susceptible or infected classes, nor move back and forth between the two classes. Thus, 


dR/dt is not coupled to the other two. 


c) 


dI/dt di BSI | pl _ iv # 
dS/dt dS —BSI' BSI BS" 
dl m 
=> =-14—. 
as + 39 
Integrating 
Lb 
T=] -1 nee 
fe i dS + ag@5: 
we get 


=> 1(S)=-S+ ie In(S) + kK where K is constant. 
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2.9 Application 9 - Creutz-Jakob Disease 


Creutz-Jakob and similar prion diseases may result from misfolded protein that 
“infects” native protein by causing it, too, to misfold, such diseases (also known as “mad 
cow disease”) lead to severe brain damage and death. Consider the schematic diagram 
below. Assume that the misfolding, like an infection, takes place when the native and 


misfolded protein come into contact. [SEK13] 


a) Propose a model for this process based on the schematic diagram. 
b) Reduce the model to a dimensionless formulation. What parameter (grouping) would 
determine whether prions are maintained or cleared from the body? 
c) Simulate the model and explore conditions under which the disease is resolved versus 


grows continually. 


ANSWER: 
a) 


0) 
Native Protein C(t) Misfolded Protein P(t) 
) ‘ 
New cell Dead cell 


Figure 2.4: Schematic Diagram of the Contact between Native and Misfiled Protein 


Again, flow outward is negative and flow inward is positive. We will assume 


the laws of mass action apply. Thus, we have 


OS = SOP ARG. 
dt 

dP 

—_— P—vP—-vP. 
‘7 oC V ay 


Let us assume that the protein population consists of only the Native and 
Misfolded proteins. Then A is the rate of new Native protein production, 7 is the rate 


of protein death, and ¢ is the rate of mass action between Native protein and Misfolded 
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protein. 


Similar to Application 7, the dimension of dC/dt and dP/dt is still [umber 
of protein cells]/|time]. The dimension of v is units of 1/time, and ¢, A, y have units of 


per protein per time. Defining x«*(t) and y*(t) as before, the dimensionless variables are 


P 
y = ave t* = pt where M = C(t) + P(t). 


Thus, 


* 


C=aMy, P=]aMs t= =. 
V 
Substituting these into a) we get 
d(y*M 
BN = ya M ~ 6ly"M)(2"M) + AWM), 


d(a*M) 


dv) = o(y*M)(a*M) — va*M — 7(2*M). 


Multiplying each term by iw we get 


dy* ee ee ee oe 
dt* ‘vp a vo 


Removing the stars, we have the dimensionless form of the equations. They are 


dy _ oM » 
dt =x ( V jay + ye? 
dx ae y 
ae ae 


The parameter (grouping) that would determine whether prion’s are maintained or 


cleared from the body is Ri = (@¢M)/v. 


c) XPP File, “diseases.ode”, to simulate the model: 
C’ = nu * P— phix C * P+ lambda* C 
P’ =—nu* P+ phi* C « P— gamma * P 
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par nu = 0.1, phi = 0.2, b = 0.05, gamma = 0.05 

init C=1,P=0 

@ xp = C, yp = P, xlo = 0, xhi = 2, ylo = —0.2, yhi = 1.2 
done 


The simulation of the model and disease conditions is: 


12 T T T T T T T T T 

og a 4 

0.6 e | 
b 

04 ae J 


Figure 2.5: XPP Simulation with Curve a) C=0.2,P=0.8 Curve b) C=0.5,P=0.5 Curve 
c) C=0.8,P=0.2 Curve d) C=1,P=0 


As you can see from curves a and 6 in the graph, when a high number of 
native protein is infected with the misfolded protein the disease will resolve itself. This 
occurs because most of the native population will already be infected. However, as seen 
with curve c, when the number of infected native protein is low the disease will grow 


continuously to infect as many native proteins as possible. 


2.10 Application 10 - Transcritical Bifurcation 


a) Redraw f(x) versus x as in Figure 1.16 with Ro = 1/2, and with Rp = 2 on the same 
set of axes. Show that Equation 2.34 has a transcritical bifurcation at 75, = 71. 


b) Use XPP Auto to redraw Figure 1.18 for —0.2 < x < 0.7 to confirm the presence of 
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a transcritical bifurcation. 

c) A student says “since the steady state x1 ceases to exist, this must be a fold bifurca- 
tion”. Explain the fallacy. [SEK13] 

ANSWER: 


a) For a transcritical bifurcation: 

1) f(x) undergoes transitions as r varies. 

2) There are two steady states. 

3) For r = 0 these graphs of f(a) merge and exchange stability. 

From the graphs seen below, it is clear that Equation 2.34 has a transcritical bifurcation 


at ss = V1. 


Figure 2.6: f(x) versus x with R, = 1/2 and with R, = 2 on the Same Set of Axes 
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Figure 2.7: Transcritical Bifurcation with —0.2 < x < 0.7 


c) The steady state x; does not cease to exist it shifts based on the value of Ro. 
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Chapter 3 


Case Study: Hepatitis C Virus 


Now that we have completed a number of exercises, we have a good under- 
standing of how a mathematical model can be created and manipulated to interpret a 
biological event. Let’s put this knowledge to use by creating a mathematical model for 
the well known Hepatitis C virus. Hepatitis C is one of six major recognizable viral 
types that cause infection of the liver. The three most common of these types are Hep- 
atitis A, Hepatitis B, and Hepatitis C. Hepatitis A virus survives in fecal matter and 
is primarily transmitted through sexual contact. The Hepatitis A virus causes only a 
short span of infection and, therefore, does not become chronic. People with Hepatitis 
A improve without treatment. The Hepatitis B virus is spread similarly to HIV, but is 
100 times more infectious because the virus can survive outside the host for many days. 
Hepatitis B can become chronic in approximately 6 percent of infected individuals and 
causes extensive damage to the liver. However, Hepatitis B symptoms are more severe 
and treatment is successful early. Hepatitis C is a bloodborne pathogen and is trans- 
mitted primarily by exposure to blood through the skin, such as through Intravenous 
Drug Use (IDU), by long-term hemodialysis, or by healthcare workers after possible 
exposure to Hepatitis C positive blood. There is no vaccine for Hepatitis C. However, 
a combination of harm reduction strategies such as the provision of new needles and 
syringes, the treatment of substance abuse, and the following of infection control guide- 
lines in healthcare are becoming a successful prevention campaign due to the growth of 


community planning groups. [AC09] 
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3.1 Hepatitis C Background 


3.1.1 Genotype 


Hepatitis C is a small, enveloped, single-stranded, positive-sense RNA virus 
that enters the cells of the liver and replicates its RNA strand in them. There are 
seven major genotypes of the Hepatiits C virus. The genotypes are divided into several 
subtypes. In the US, about 70 percent of the cases are cause by genotype 1 and 20 
percent by genotype 2. The incubation period for the Hepatitis C virus ranges from 
14-168 days with an average of 28-48 days. A combination of tests is required to diag- 
nose Hepatitis C. The initial test is the Hepatitis antibody enzyme immunoassay test 
which indicates either past or present infection. If the Hepatitis antibody enzyme test 
is positive, a Polymerase Chain Reaction (PCR) test is given. This test detects the 
presence of the Hepatitis C virus in the blood and, thus, is used to diagnose chronic 


Hepatitis C infection. [HSAH98] 


3.1.2. Hepatitis C Virus Progression 


The initial symptoms of Hepatitis C are typically misinterpreted and often 
disappear after a few weeks. This initial stage of the condition is known as Acute 
Hepatitis C. During the acute infection period, the body’s immune system is at work 
trying to kill the Hepatitis C infected cells. In 15 to 25 percent of the cases, the body 
is successful and the individuals do not suffer further infections. They are considered 
“cured” or “dormant” if they remain infection free for 6 months. For the other 75-85 
percent, the infections continue and the accumulation of scar tissue resulting from these 
infections eventually prevents the liver from functioning properly. In approximately 15 
percent of these cases, a person develops cirrhosis, necrosis, and, then, cancer. Five 
percent die. This stage of the Hepatitis C progression is known as Chronic Hepatitis 
C and is defined as infection with the Hepatitis C virus recurring for more than six 
months based on the presence of the Hepatitis C single-stranded RNA in the blood. 
Chronic infections are commonly asymptomatic during the first two decades and, thus, 
are consistently discovered following the investigation of elevated liver enzyme levels or 
during a routine screening of high-risk individuals. There is no cure for Hepatitis C, but 


a proper medication regime can stop the virus from replicating itself. [HSAH98] [AC09] 
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3.1.3 Treatment 


The primary goal of treatment is to achieve a Sustained Viral Response (SVR), 
which is defined as undetectable Hepatitis C virus in the blood 6 months after the end 
of treatment. Prior to 2011, treatments consisted of a combination of pegylated in- 
terfon alpha and ribavirin for a period of 24-48 weeks depending on the Hepatitis C 
virus genotype. This combination provided SVR “cure” rates between 70-80 percent for 
patients with genotype 2 and 3 and 45-70 percent for patients with genotype 1 and 4. 
After 2011, treatments of a combination of sofosbuvir with ribavirin and interfon have 
proved to be approximately 90 percent effective in patients with genotype 1,4,5, and 6 
and treatments of sofosbuvir with only ribavirin have proved to be approximately 70-95 


percent effective in patients with genotype 2 and 3. [HSAH98] 


An estimated 4 million Americans (1.8 percent of all Americans) have been 
infected with the Hepatitis C virus. Based upon national data, in 2012, an estimated 
600,000 Californians are currently infected with Hepatitis C and 5,000 Californians are 
newly infected each year. [HSAH98] 


3.2 Hepatitis C Virus Model, Equations, and Analysis 


3.2.1 Block Diagram of Hepatitis C Virus Progression 


The schematic (Block) diagram for the Hepatitis C virus progression is as follows 


[HVCACC11]: 


Figure 3.1: Schematic Diagram - Hepatitis C Virus Progression. 


We define: 

S(t) = the number of Susceptible people 
A(t) = the number of Acute infected people 

C(t) = the number of Chronic infected people 

D(t) = the number of Dormant infected people 

L(t) = the number of infected people with Liver disease 


N(t) = the number of infected people with Necrosis of the Liver 


( = the rate increase of infected people due to IVDA (intravenous drug abuse) 
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6 = the rate increase of chronic infected people due the bodies inability to kill the active 


hepatitis C cells 


6 = the rate increase of acute infected people with sustained viral response (6months) 


p = the rate increase of hepatitis C dormant people due to medication, diet, and alcohol 
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abstinence 


a = the rate increase of people with unsustained response (6 months) 
pu = the rate increase of extended chronic hepatitis C infected people 


w = the rate increase of people with extended liver disease 


Then the equations are: 


dS 
ap Pe: 
a = BSA-—O0A-— dA, 
° =0A+ pC — aD, 
“ = d6A+aD — pC — uC, 
“ = pC —wl, 
=wl. 
In this case, the total conserved population is au = as 4 dA | ie ae + a ix = 


BSA+PBSA-O0A-—6A+60A+ pC —aD+6A+aD— pC —pC+pC—wlh+wLl=0. 


3.2.2 Simplified Model of Hepatitis C Virus 


To simplify the model into a set of ordinary differential equations that we have 
discussed so far, we will consider only the relationship between the Dormant infected 
individuals and the Chronic infected individuals. Again, once an individual is infected 
with the Hepatitis C virus, the virus will always be present. Therefore, the infected pop- 
ulation, I(t), will either have chronic Hepatitis C infections, C(t), for the next couple 
of decades or they will be considered “cured” (virus is dormant) by a Sustained Viral 
Response (SVR) as described above. Thus, we have the following schematic (block) 


diagram and differential equations: 
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c(t) [) «~Dit) 


Figure 3.2: Schematic Diagram of Chronic and Dormant Hepititis C Infected Individuals 


ee aD 
oe ee 
dD 
a ee 
with I(t) = D(t) + C(t) thus, 
dl dD dd 
= =—4 D+aD = 0, 
to a RE 


so conservation is preserved and I is constant. 


Therefore, substituting [(t) = D(t) + C(t) into the equations gives us 


- =a(I—C) — pC, (1) 
7 =plI-D)-aD, (2) 


dC and we have dimensions of (number of people) / (time) and p and a have units of 


dt 
1 / (time) (days~'). 


Defining x*(t) as the fraction of the number of infected in the Chronic stage 


and y*(t) as the fraction of the number of infected in the Dormant stage. 


We define dimensionless variables 


x T’ ts, er 


Since, C(t) + D(t) = I = constant, it follows that C/I + D/I=a*+y* =I/I=1. 
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Substituting, C= 2*I, D=y*I, t =t*/p into (1) and (2) we have 


da") owt — oat 
“Co a le 
d(y*I) _ : i 
d/o) PE -¥ I) -ay*l. 


Dividing both equations by J and p, we get 


= —(l-2*)- 2’, 
a8 a ) 
dy* a Qa 
= 1 * ok = 1 _ * 1 =) 
AG y y eil+—} 
Dropping the stars, we have 
ae (l-—a)-2 
dtp , 
dy ay 
=1 14 
Ai y( Fe 
Now, we will find the value of x at steady-state dx/dt = 0. Solving for x,,, we 
obtain 
d 
a= (i x)—-x=0, 
dt =p 
=> 2-"y-25 0, 
pp 
=> —-—2(—+1)=0, 
a 
— >! pP = 
Bek. OAD 
= ray 
= +4 bf = 
8s atp 
To convert back to unit-carrying variables, we multiply by J, so that C = xI, 
then 
eg 
a+ p 


Thus, if = a/(a +p) is reached, the virus returns to an active state and the 


patient becomes chronic. Notice, this steady state is feasible only if « > 0 which means 
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that ; > 0. Therefore, : must be positive for the virus to progress to the chronic state. 


Figure 3.3 below shows the sketch of 1 = %—2(5 +1) = f(x). As seen, f(z) is 


a straight line with arrows showing the direction of change of x(t) from various initial 


values. 


Figure 3.3: Sketch of f(z) vs. x a= .64,p = .36. 


The model simulations shown in Figure 3.4, 3.5, and 3.6 below are plots of x(t) 
versus t for various initial conditions. The first and second simulations show the fraction 
of the infected population that are in the chronic stage of Hepatitis C before 2011 with 
and without the influence of diet and alcohol abstinence. In this case, the available medi- 
cations consisted of only the combination of Pegylated Interfon Alpha and Ribavirin. As 
you can see, there is a significant decrease of chronic infected individuals when diet and 
alcohol abstinence are maintained. The third simulation shows the fraction of the in- 
fected population that are in the chronic stage of Hepatitis C currently. In this case, the 
available medications consist of a combination of Sofosbuvir with Ribavirin and Interfon. 
The values of a and p used in these three simulations were calculated from the percent- 
age of infected individuals that will become chronic, the percentage of chronic infected 
individuals with Hepatitis C genotype 1 and 2, and the percentage of chronic infected in- 
dividuals that respond to medication with and without proper diet and abstinence from 


alcohol. Therefore, Simulation 1, before 2011 and without diet and alcohol abstinence, 
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p = .70(.80).45 + .20(.80).70 = .36 and a = .64. Simulation 2, before 2011 and with diet 
and alcohol abstinence, p = .70(.80).50 + .20(.80).80 = .41 and a = .59. Simulation 3 
after 2011 and without diet and alcohol abstinence, p = .70(.80).90 + .20(.80).70 = .62 
and a@ = .38. The difference between the first two fractions of infected individuals in 
Simulation 1 and Simulation 2 can be shown to fluctuate depending on the number of 
chronic infected individuals who follow the guidelines of diet and alcohol abstinence. 
Because of this, over the last decade, healthcare professionals and community planning 
groups are increasing the availability and amount of education programs given to chronic 
infected individuals on the importance of diet and alcohol abstinence. The difference 
between the two fractions of infected individuals in Simulation 1 and Simulation 8 is di- 
rectly proportional to the advancement in medications available to Hepatitis C infected 
individuals. As you can see, the advanced treatment has greatly decreased the fraction 


of the population with known chronic Hepatitis C. 


Figure 3.4: Simulation 1 (Infection rates prior to 2011 without Diet and AA) a = 
64, p = .36. 
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Figure 3.5: Simulation 2 (Infection rates prior to 2011 with Diet and AA) a = .59,p = 
Al. 
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Figure 3.6: Simulation 3 (Infection rates after 2011 without Diet and AA) a = .38,p = 
62. 
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3.2.3 Future Studies-The Next Step 


Remember, we kept the model and equations limited to a simple, single or- 
dinary differential equation for the purpose of remaining within the boundaries of the 
material reviewed in this paper. The next step would be to apply one more variable to 
the model and increase the complexity to a system of two first-order differential equa- 


tions. In doing so, we would use Phase Plan methods to study these systems, generally, 


written as: 
dx 
dt i Fa; Y), 
dy 
aa g(x,y). 


To show an example of a system of two first-order differential equations, let us 


add a variable to the simple Hepatitis C schematic (block) diagram: 


Figure 3.7: Schematic Diagram - Susceptible, Dormant, and Infected Hepitits C 
Individuals 


Let I(t) will represent the combination of A(t) and D(t), the total population 


of infected individuals that are not chronic. The differential equations for this model 


are: 
dS 
— = -BSI - BSC, 
dl 
7G 7 OSI—al + pC + BSC, 
dC 


Total population is conserved: 


dN aS adI dC 
= = I I I4 Lol =0. 
‘7 . + di BS BSC + BS a pC+6SC+a pC =0 
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Therefore, N is constant and S = N—I—C. Substituting this equation for all S in the 


differential equations, we get: 


“= 6(N —I-C)I-B(N-I-C)C, 

dl 

GH BIN -1-C) ~ al + pC + BIN -1-C)C, 
dC 
Be ap ape! 
et 


which is a system of equations with two unknowns, namely, J and C. 

From here, solutions to these differential equations are found by graphing an 
xy plane (phase plane) which shows all trajectories of the solutions through every point. 
Steady-states are then found by finding the intersection of the x and y nullcline curves 


da _ dy _ 
of G = 0 and it = 
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